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We describe a maximum entropy approach for computing volumes and counting in- 
teger points in polyhedra. To estimate the number of points from a particular set 
X C M n in a polyhedron P C M™, by solving a certain entropy maximization prob- 
lem, we construct a probability distribution on the set X such that a) the probability 
mass function is constant on the set PC\X and b) the expectation of the distribution 
lies in P. This allows us to apply Central Limit Theorem type arguments to deduce 
computationally efficient approximations for the number of integer points, volumes, 
and the number of 0-1 vectors in the polytope. As an application, we obtain as- 
ymptotic formulas for volumes of multi-index transportation polytopes and for the 
£SJ ■ number of multi-way contingency tables. 
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i/-) ■ 1. Introduction 



In this paper, we address the problems of computing the volume and counting 
the number of integer points in a given polytope. These problems have a long his- 
tory, see for example, surveys [GK94], [DL05] and [Ve05], and, generally speaking, 
are computationally hard. We describe a maximum entropy approach which, in a 
number of non-trivial cases, allows one to obtain good quality approximations by 
solving certain specially constructed convex optimization problems on polytopes. 
Those optimization problems can be solved quite efficiently, in theory and in prac- 
tice, by interior point methods, see [NN94]. 

The essence of our approach is as follows: given a discrete set S C M. n of interest, 
such as the set of all no n- negative integer points or the set {0, l} n of all 0-1 
points, and an affine subspace A C M. 71 we want to compute or estimate the number 
1 5 D A\ of points in A. For that, we construct a probability measure \i on S with 
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the property that the probability mass function is constant on the set A n S and 
the expectation of ji lies in A. These two properties allow us to apply Local Central 
Limit Theorem type arguments to estimate \S n A\. The measure \i turns out to 
be the measure of the largest entropy on S with the expectation in A, so that 
constructing [i reduces to solving a convex optimization problem. We also consider 
a continuous version of the problem, where S is the non-negative orthant IR™ and 
our goal is to estimate the volume of the set S n A. 

Our approach is similar in spirit to that of E.T. Jaynes [Ja57] (see also [Go63]), 
who, motivated by problems of statistical mechanics, formulated a general princi- 
ple of estimating the average value of a functional g with respect to an unknown 
probability distribution on a discrete set S of states provided the average values of 
some other functionals f\ , . . . , f r on S are given. He suggested to estimate g by 
its expectation with respect to the maximum entropy probability distribution on S 
such that the expectations of fi have prescribed values. Our situation fits this gen- 
eral framework when, for example, S is the set Z" of non-negative integer vectors, 
fi are the equations defining an affine subspace A, functional g is some quantity of 
interest, while the unknown probability distribution on S is the counting measure 
on S H A (in interesting cases, the set S H A is complicated enough so that we may 
justifiably think of the counting measure on S H A as of an unknown measure). 

(1.1) Definitions and notation. In what follows, M n is Euclidean space with 
the standard integer lattice Z n C M n . A polyhedron P C W 1 is defined as the set 
of solutions x = (£1, . . . , £ n ) to a vector equation 

(1.1.1) iiai + . . . + i n a n = 6, 

where a\, . . . , a n ; b G ~R d are <i-dimensional vectors for d < n, and inequalities 

(1.1.2) > o. 

We assume that vectors ai,... , a n span R d , in which case the affine subspace 
defined by (1.1.1) has dimension n — d. We also assume that P has a non-empty 
interior, that is, contains a point x = , £ n ), where inequalities (1.1.2) are 

strict. One of our goals is to compute the (n — d) -dimensional volume volP of P 
with respect to the Lebesgue measure in the affine subspace (1.1.1) induced from 
M n . More generally, our approach allows us to estimate the exponential integral 

J e 1 ^ dx, 

where I : R n — >■ R is a linear function. We note that the integral may be well 
defined even if P is unbounded. Often, we use a shorthand Ax = 6, x > for 
(1.1.1)— (1.1.2), where A = [ai, . . . , a n ] is the matrix with the columns a\, . . . ,a n 
and x is thought of as a column vector x — . . . , £ n ] . 
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We are also interested in the number |PflZ n | of integer points in P. In this case, 
we assume that vectors a\, . . . ,a n and b are integer, that is, a\, . . . , a n ; b G 7L d . 
The number \P n Z n | as a function of vector b in (1.1.1) is known as the vector 
partition function associated with vectors a\, . . . ,a n , see for example, [BV97]. More 
generally, our approach allows us to estimate the exponential sum 

mePnz™ 

where £ : R n — > 1 is a linear function. Again, the sum may converge even if 
polyhedron P is unbounded. 

Finally, we consider a version of the integer point counting problem where we 
are interested in 0-1 vectors only. Namely, let {0, l} n be the set (Boolean cube) of 
all vectors in R n with the coordinates and 1. We estimate |Pn{0, l} n | and, more 
generally, the sum 

m£Pn{0,l}" 

(1.2) The maximum entropy approach. Let us consider the integer counting 
problem first. One of the most straightforward approaches to computing \P D Z n | 
approximately is via the Monte Carlo method. As in Section 1.1, we think of P as 
defined by a system Ax = 6, x > 0. We place P in a sufficiently large axis-parallel 
integer box B in the non-negative orthant R™ of R n , sample integer points from 
B independently at random and count what proportion of points lands in P. It is 
well understood that the method is very inefficient if P occupies a small fraction of 
£?, in which case the sampled points will not land in P unless we use great many 
samples. Let X be a random vector distributed uniformly on the set of integer 
points in box B. One can try to circumvent sampling entirely by considering the 
random vector Y = AX and interpreting the number of integer points in P in terms 
of the probability mass function of Y at b. One can hope then, in the spirit of the 
Central Limit Theorem, that since the coordinates of Y are linear combinations 
of independent coordinates xi, . . . ,x n of X, the distribution of Y is somewhat 
close to the Gaussian and hence the probability mass function of Y at b can be 
approximated by the Gaussian density. The problem with this approach is that, 
generally speaking, the expectation E Y will be very far from the target vector b, so 
one tries to apply the Local Central Limit Theorem on the tail of the distribution, 
which is precisely where it is not applicable. 

We propose a simple remedy to this naive Monte Carlo approach. Namely, by 
solving a convex optimization problem on P, we construct a multivariate geometric 
random variable X such that 

(1.2.1) The probability mass function of X is constant on the set P fl Z n of 
integer points in P; 



(1.2.2) We have EX e P, or, equivalently, EY = b for Y = AX. 
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Condition (1.2.1) allows us to express the number \P fl Z n | of integer points in 
P in terms of the probability mass function of Y, while condition (1.2.2) allows 
us to prove the Local Central Limit Theorem for Y in a variety of situations. We 
have X = (xi, . . . ,x n ) where Xj are independent geometric random variables with 
expectations Q such that z = (Ci, • • • , Cn) is the unique point maximizing the value 
of the strictly concave function, the entropy of X, 

n 

0(*) = £(& + i)Mfc + i)-&K*) 

on P, see Theorem 3.1 for the precise statement. 

Similarly, to estimate the number of 0-1 vectors in P, we construct a multivariate 
Bernoulli random variable X, such that (1.2.2) holds while (1.2.1) is replaced by 

(1.2.3) The probability mass function of X is constant on the set P fl {0, l} n of 
0-1 vectors in P. 

In this case, X = (xi, . . . , x n ), where Xj are independent Bernoulli random vari- 
ables with expectations Q such that z = (Ci, • • • , Cn) is the unique point maximizing 
the value of the strictly concave function, the entropy of X, 

Kx) = J2 (& ln } 3 + c 1 - ^ ln rr^) 

on the truncated polytope 

Pn{0<e,<l: for j = l,...,n}, 

see Theorem 3.3 for the precise statement. 

Finally, to approximate the volume of P, we construct a multivariate exponential 
random variable X such that (1.2.2) holds and (1.2.1) is naturally replaced by 

(1.2.4) The density of X is constant on P. 

Condition (1.2.4) allows us to express the volume of P in terms of the density 
of Y = AX at Y = 6, while (1.2.2) allows us to establish a Local Central Limit 
Theorem for Y in a number of cases. In this case, each coordinate Xj is sam- 
pled independently from the exponential distribution with expectation Q such that 
z = (Ci, • • • , Cn) is the unique point maximizing the value of the strictly concave 
function, the entropy of X, 



n 

f( x ) =n + 5^1n£ J - 



on P, see Theorem 3.6 for the precise statement. In optimization, the point z is 
known as the analytic center of P and it played a central role in the development 
of interior point methods, see [Re88]. 

These three examples (counting integer points, counting 0-1 vectors, and com- 
puting volumes) are important particular cases of a general approach to counting 
through the solution to an entropy maximization problem (cf. Theorem 3.5) with 
the subsequent asymptotic analysis of multivariate integrals needed to establish the 
Local Central Limit Theorem type results. 

2. Main results 

(2.1) Gaussian approximation for volume. Let F C 1" be a polytope, defined 
by a system Ax = 6, x > 0, where A is an d x n matrix with the columns ai, . . . , a n . 
We assume that rank A = d < n. We find the point z = (Ci, • • • , Cn) maximizing 

n 

f{x) = n + ^ln^-, x= (fi,... ,£n) 
i=i 

on P. Let B be the dxn matrix with the columns Ci a i? ■ ■ ■ ■> Cn a n- We approximate 
the volume of P by the Gaussian formula 

1 fdetAA T \ 1/2 f(z) 
< 2 - L1 > VOlP ^(2^(detWj 

We consider the standard scalar product (•, •) and the corresponding Euclidean 
norm || • || in M d . 

We prove the following main result. 

(2.2) Theorem. Let us consider a quadratic form q : M. d — > R defined by 

1 n 

i=i 

Suppose that for some A > we have 

q(t) > X\\tf for all t G R d 

and that for some 9 > we have 

CjhjW <0 for j = 1, . . . ,n. 

Then there exists an absolute constant 7 such that the following holds: 
let < e < 1/2 be a number and suppose that 
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A > 7 # 2 e - 2 (Wln^ ln(^) 
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Then the number 

1 ( detAA- \ 1/2 f(z) 
(2n) d / 2 Vdet BB T ) 
approximates volP within relative error e. 

Let us consider the columns ai, . . . , a n of A as vectors from Euclidean space W 1 
endowed with the standard scalar product (•,•). The quadratic form q defines the 
moment of inertia of the set of vectors {Ci«i, • • • , Cn«n} 5 see, for example, [Ba97]. 
By requiring that the smallest eigenvalue of q is sufficiently large compared to the 
lengths of the vectors C,j a ji we require that the set is sufficiently "round". For a 
sufficiently generic (random) set of n vectors, we will have q(t) roughly proportional 
to ||t|| 2 and hence A will be of the order of nd~ l maxj = i r .. >n C|||°il| 2 - 

We prove Theorem 2.2 in Section 6. 

In Section 4, we apply Theorem 2.2 to approximate the volume of a multi-index 
transportation polytope, see, for example, [Y+84], that is, the polytope P of v- 
dimensional k\ x . . . x k v arrays of non-negative numbers with 1 < ji < ki 

for i = 1, . . . , v with prescribed sums along coordinate hyperplanes ji = j. We show 
that Theorem 2.2 implies that asymptotically the volume of P is given by a Gaussian 
formula (2.1.1) as long as v > 5. We suspect that the Gaussian approximation holds 
as long as v > 3, but the proof would require some additional considerations beyond 
those of Theorem 2.2. In particular, for v > 5 we obtain the asymptotic formula 
for the volume of the polytope of polystochastic tensors, see [Gr92]. 

For v = 2 polytope P is the usual transportation polytope. Interestingly, its 
volume is not given by the Gaussian formula, cf. [CM07b] . 

In [Ba09] , a much cruder asymptotic formula of vol P in terms of e^ z > was proved 
under much weaker assumptions. 

(2.3) Gaussian approximation for the number of integer points. For a 

polytope P, defined by a system Ax = 6, x > 0, we find the point z = (Ci, • • • , Cn) 
maximizing 

n 

9(x) = (fo + !) ln + !) " . * = (fr, . . . , £ n ) 
on P. Assuming that a i , . . . , a n G Z are the columns of A, we define B as the 

1/2 

d x n matrix whose j-th column is (£j + Cj) a,j for j = 1, ... ,n. 

We assume that A is an integer d x n matrix of rank d < n. Let A = A (Z n ) 
be image of the standard lattice, A C Z d . We approximate the number of integer 
points in P by the Gaussian formula 

e fO) det A 

(2.3.1) \PDZ nl ' 



(27r) d / 2 (detSS T )V2- 



In this paper, we consider the simplest case of A = Z d , which is equivalent to the 
greatest common divisor of the d x d minors of A being equal to 1. 
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Together with the Euclidean norm || • || in M d , we consider the i 1 and £°° norms: 

d 

Plli = / l r j| an d \\t\\oo = max \ti\ where t = (ri, . . . , r^) . 

' i=l,... ,d 



Clearly, we have 



l^lli > \\t\\ > Plloo for all t £ 



Compared to the case of volume estimates (Sections 2.1-2.2), we acquire an additive 
error which is governed by the arithmetic of the problem. 

Let ei, . . . , ea be the standard basis of Z d . We prove the following main result. 

(2.4) Theorem. Let us consider a quadratic form q : M. d — > R defined by 

1 n 

For i = 1, . . . ,d let us choose a non-empty finite set Y\ C Z n swc/i t/iat Ay = /or 
a// y G Yi and let us define a quadratic form ipi : R n — ► R &?/ 



Suppose that for some X> we have 

q(t)>X\\tf for all teR d , 
that for some p > we have 

ipi{x) < p\\x\\ 2 for all x G R n and i = l,...,d, 
that for some 9 > 1 we have 



and taat 

0(1 + 0) > « /or J = 1,... ,n 

and some a > 0. 

Then, for some absolute constant 7 > and /or any < e < l/2 ; as Zona as 



A>7e- 2 # 2 (d + ln^ 2 ln(^), 
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we have 



Pf]Z n \ = 



K 



+ A 



(27r) d / 2 (det BB T ) 



1/2 



where 



l-e< K<l + e 



and 




) 



— m 



1 



for m = 



16ttV 2 



While the condition on the smallest eigenvalue of quadratic form q is very similar 
to that of Theorem 2.2 and is linked to the metric properties of P, the appearance 
of quadratic forms fa is explained by the arithmetic features of P. Let us choose 
1 < % < d and let us consider the affine subspace Ai of the points x G M n such 
that Ax = ej. Let Aj = Ai fl Z n be the point lattice in Ai. We would like to 
choose a set lj C Aj in such a way that the maximum eigenvalue pi of the form fa, 
which defines the moment of inertia of lj, see [Ba97], becomes as small as possible, 
Pi <C 1, so that the additive error term A becomes negligibly small compared to 

1/2 

the Gaussian term (27r) _d / 2 (detSS T ) . For that, we would like the set lj to 
consist of short vectors and to look reasonably round. Let us consider the ball 
B r = {x G K n : ||x|| < r} of radius r and choose Fj = B r fl Aj. If the lattice points 
Yi are sufficiently regular in B r fl Ai then the moment of inertia of Y± is roughly the 
moment of inertia of the section B r fl Ai, from which it follows that the maximum 
eigenvalue of fa is about r 2 / dim^4j = r 2 /(n — d). Roughly, we get 



where r is the smallest radius of the ball B r such that the lattice points B r fl Aj 
are distributed regularly in every section B r fl Ai for i = 1, . . . ,d. 
We prove Theorem 2.4 in Section 8. 

In Section 5, we apply Theorem 2.4 to approximate the number of 1-margin 
multi-way contingency tables, see for example, [Go63] and [DO04], that is, v- 
dimensional k\ x . . . x k v arrays of non- negative integers (£j lm ..j v ) with 1 < jj < k{ 
for % = 1, . . . , v with prescribed sums along coordinate hyperplanes ji = j. We 
show that Theorem 2.4 implies that asymptotically the number of such arrays is 
given by a Gaussian formula (2.3.1) as long as v > 5. We suspect that the Gaussian 
approximation holds as long as v > 3, but the proof would require some additional 
considerations beyond those of Theorem 2.4. 

In [Ba09], a much cruder asymptotic formula with the main term e 9 ^ in the 
logarithmic order is shown to hold for the number of integer points in flow polytopes 
(a class of polytopes extending transportation polytopes for v = 2). At our request, 
A. Yong [Yo08] computed a number of examples. Here is one of them, originating in 
[DE85] and then often used as a benchmark for various computational approaches: 



P 



(n-dY 
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we want to estimate the number of 4 x 4 non-negative integer matrices with row 
sums 220, 215, 93 and 64 and column sums 108, 286, 71 and 127. The exact number 
of such matrices is 1225914276768514 m 1.23 x 10 15 . Framing the problem as the 
problem of counting integer points in a polytope in the most straightforward way, we 
obtain an over-determined system Ax = b (note that the row and column sums of a 
matrix are not independent). Throwing away one constraint and applying formula 
(2.3.1), we obtain 1.30 x 10 15 , which overestimates the true number by about 6%. 
The precision is not bad, given that we are applying the Gaussian approximation to 
the probability mass-function of the sum of 16 independent random 7-dimensional 
integer vectors. 

(2.5) Gaussian approximation for the number of 0-1 points. For a polytope 
P defined by a system Ax = 6, < x < 1 (shorthand for < < 1 for x = 
(£i 5 • • • 7 £n)), we find the point z = . . . , ( n ) maximizing 

n i i 

h{x) = J2 ln r + (1 " & ln r^r) ' x = • • • ' 60 

3 = 1 ^ ^ J 

on P. Assuming that A is an integer matrix of rank d < n with the columns 
a\, . . . ,a n G Z d , we compute the d x n matrix B whose j-th column is 

1/2 

(Cj — Cj) a j- We approximate the number of 0-1 vectors in P by the Gauss- 
ian formula 

e h ( z ) Het A 

(2.5.1) |Pn{o,i} n | 



(2 7 r) d / 2 (detSS T ) 1 /2' 



where A = A (Z n ). Again, we consider the simplest case of A = Z d . We prove the 
following main result. 

(2.6) Theorem. Let us consider a quadratic form q : M. d — > R defined by 

1 n 

?(*) = 2Efo-Ci 2 ) <^> 2 - 

For i = 1, . . . ,d let us choose a non-empty finite set Yi C Z n such that Ay = for 
all y EYi and let us define a quadratic form ipi : R n — >■ R by 



\ i\ r-V 



yeYi 



Suppose that for some X > we have 



q(t)>X\\t\\ 2 for all t G M. d , 
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that for some p > we have 



ipi(x) < p\\x\\ for all x £ W 1 and i = l, 
that for some 9 > 1 we have 



hjWi < flyOCi-Cj) f° r j = i,...,n 

and that for some < « < 1/4 roe have 

0(1-0) > « /or j = 1, ... ,n. 
Then, for some absolute constant 7 > and for any < e < 1/2, as long as 

A> 7 e" 2 # 2 (d + ln^ln^), 

we have 

\p n {0, l} n | = e^ 2 ) f + A 

V(27r) d / 2 (det BB T ) 1/2 

where 

1-e < k <l + e 

a 



|A| < exp 



80^ 2 pJ ' 



We note that in [Ba08] a much cruder asymptotic formula with the main term 
e h(z) j n logarithmic order is shown to hold for the number of 0-1 vectors in flow 
polytopes. 

We prove Theorem 2.6 in Section 7. 

In Section 5, we apply Theorem 2.6 to approximate the number of binary 1- 
margin multi-way contingency tables, see for example, [Go63] and [DO04], that is, 
z/-dimensional k± x . . . x h v arrays (^j 1 ...j v ) of 0's and l's with 1 < ji < h\ for 
% = 1, . . . , v with prescribed sums along coordinate hyperplanes ji = j. Alterna- 
tively, the number of such arrays is the number of z/-partite uniform hypergraphs 
with prescribed degrees of all vertices. We show that Theorem 2.6 implies that 
asymptotically the number of such arrays is given by the Gaussian formula (2.5.1) 
as long as v > 5. We suspect that the Gaussian approximation holds as long as 
v > 3, but the proof would require some additional considerations beyond those of 
Theorem 2.6. 
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3. Maximum entropy 



We start with the problem of integer point counting. 

Let us fix positive numbers p and q such that p + q = 1. We recall that a discrete 
random variable x has geometric distribution if 

Pr [x = k] = pq k for k = 0, 1, . . . . 

For the expectation and variance of x we have 

E x = — and var x = 

p p z 

respectively. Conversely, if E x = Q for some ( > then 

1 C 

p = Y^' q = Y+~c and vara; = £ + £ 2 - 

Our first main result is as follows. 

(3.1) Theorem. Let P C M. n be the intersection of an affine subspace in R n and 
the non-negative orthant R". Suppose that P is bounded and has a non-empty 
interior, that is contains a point y = (771 , . . . , rj n ) where rjj > for j = 1, . . . , n. 
Then the strictly concave function 

n 

^) = 5]((0 + l)ln(0 + l)-01ne,) for x = (£1, . . . , £ n ) 

attains its maximum value on P at a unique point z = (£1, . . . , £ n ) sitc/i £/iat (j > 
for j = 1, ... ,n. 

Suppose now that Xj are independent geometric random variables with expecta- 
tions Q for j = 1, . . . , n. Let X = (xi, . . . , x n ). Then the probability mass function 
of X is constant on P fl Z n and equal to e~ 9 ^ at every x G P fl Z n . In particular, 

\PnZ n \ = e 9{z) Pr {X e P). 



Proof. It is straightforward to check that g is strictly concave on the non-negative 
orthant R™, so it attains its maximum on P at a unique point z = (Ci>- • • ,Cn)- 
Let us show that Q > 0. Since P has a non-empty interior, there is a point 
y = (7/1, . . . , 7/ n ) with rjj > for j = 1, . . . , n. We note that 

— y = In ^- — 

which is finite for £j > and equals +00 for = (we consider the right derivative 
in this case). Therefore, if Q = for some j then ^((1 — e)z + ey) > g(z) for all 
sufficiently small e > 0, which is a contradiction. 
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Suppose that the affine hull of P is defined by a system of linear equations 

n 

Y a 'i^' = & for * = 1 » • • • ' d - 

3 = 1 

Since z is an interior maximum point, the gradient of g at z is orthogonal to the 
affine hull of P, so we have 



(i + C' \ d 



and some Ai, . . . , A^. Therefore, for any x e P, a; = (£i, . . . , £ n ), we have 

n / ^ -)- ^ ■ \ n d d 

Y & ln ( -t-^ ) = Y A '^ '/ = A *&> 

j = l ^ ^ J ' j = l i=l i=l 



or, equivalent ly, 

(3.1.1) u(^y =^[t, x ^\- 



Substituting £j = ^ for j = 1, . . . , n, we obtain 

r d 

1 = ex p \ Y 

From (3.1.1) and (3.1.2), we deduce 



(3.1.2) n( 1 ^) G =exp|x:A l A 



.7 = 1 ^ i 



+ 0) 



i+0 



The last identity states that the probability mass function of X is equal to e 

for every integer point x G P. □ 

One can observe that the random variable X of Theorem 3.1 has the maximum 
entropy distribution among all distributions on Z™ subject to the constraint ElG 
P. 
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Theorem 3.1 admits the following straightforward extension. Let I : R n — > R 
be a linear function, 

= 71C1 + • • • + 7nCn where x = . . . , f„) . 

Let P C R n be a polyhedron as in Theorem 3.1, although not necessarily bounded, 
and suppose that I is bounded on P from above and attains its maximum on P on 
a bounded face of P (it is not hard to see that this condition is sufficient for the 
series J2 x ePnz n exp{£(x)} to converge). Then the strictly concave function 

n 

9iW = E (& + !) ln & + !) - & ln & + litj) 

attains its maximum on P at a unique point z = (0, • • • , Cn)> where £j > for j = 
1, . . . , n. Suppose now that X = (xi, . . . , x n ) is the vector of independent geometric 
random variables such that Exj = Q for j = 1, . . . , n. Then the probability mass 
function of X at a point x E P fl Z n is equal to exp {— ^(z) + £(x)}. In particular, 

^ exp{£(x)} = exp{g e (z)}Pr{XeP}. 

The proof is a straightforward modification of that of Theorem 3.1. 

(3.2) The Gaussian heuristic for the number of integer points. Below we 
provide an informal justification for the Gaussian approximation formula (2.3.1). 

Let P be a polytope and let X be a random vector as in Theorem 3.1. Suppose 
that P is defined by a system Ax = 6, x > 0, where A = (oiy ) is a d x n matrix of 
rank d < n. Let Y = AX, so Y = (yi, . . . , j/d), where 

n 

yi = a^Xj for z = 1, . . . , d. 

3 = 1 

By Theorem 3.1, 

|PHZ n | = e 9{z)l PY {Y = b) 

and 

EF = Az = 6. 

Moreover, the covariance matrix Q = (g^ ) of Y is computed as follows: 

n n 

q i:j = cov (y i5 ?/j) = ^ «ifcttjfcvar x fc = a^a^ (Cfc + Cjfc) • 
k=i fc=i 

We would like to approximate the discrete random variable Y - by the Gaussian 
random variable Y* with the same expectation b and covariance matrix Q. We 
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assume now that A is an integer matrix and let A = {Ax : x G Z n |. Hence 
A C Z d is a (i-dimensional lattice. Let IT C lR d be a fundamental domain of A, so 
vol IT = det A. For example, we can choose IT to be the set of points in R d that are 
closer to the origin than to any other point in A. Then we can write 

|PnZ n | = e 9{z) Pr{Y G b + U}. 

Assuming that the probability density of Y* does not vary much on b + IT and that 
the probability mass function of Y at Y = b is well approximated by the integral 
of the density of Y* over b + IT, we obtain (2.3.1). 

Next, we consider the problem of counting 0-1 vectors. 

Let p and q be positive numbers such that p + q = 1. We recall that a discrete 
random variable x has Bernoulli distribution if 

Pr{x = 0}=p and Pr {x = 1} = q. 

We have 

E x = q and var x = qp. 
Conversely, if E x = ( for some < £ < 1 then 

p=l-C, q = C and varx = (-( 2 . 
Our second main result is as follows. 

(3.3) Theorem. Let P C W 1 be the intersection of an affine subspace in W 1 
and the unit cube {0 < £j < 1 : j = 1,... , n}. Suppose that P has a non- 
empty interior, that is, contains a point y = (771,... , rj n ) where < rjj < 1 for 
j = 1, . . . , n. Then the strictly concave function 

n ( 1 1 \ 

M^) = E^ln- + (l-0)ln T -^-J for x = (£1, . . . ,£„) 

attains its maximum value on P at a unique point z = (£1, . . . , £ n ) swc/i t/iat < 
< 1 /or j = 1, . . . ,n. 

Suppose now that Xj are independent Bernoulli random variables with expecta- 
tions Q for j = 1, . . . , n. Let X = (xi, . . . , x n ). Then the probability mass function 
of X is constant on P fl {0, l} n and equal to e~ h ^ for every xG?n{0, l} n . In 
particular, 

|Pn{0,l} n | = e h{z) Pr{X e P}. 

□ 

One can observe that X has the maximum entropy distribution among all distri- 
butions on {0, l} n subject to the constraint EIgF. The proof is very similar to 
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that of Theorem 3.1. Besides, Theorem 3.3 follows from a more general Theorem 
3.5 below. 

Again, there is a straightforward extension for exponential sums. For a linear 
function £ : R n — > R, 



£(x) = 71& + . . . + 7 n £n where x = . . . , £ n ) , 



we introduce 



n / 1 1 \ 

M*) = 5] ln f. + i 1 - ln + TjOJ • 

Then the maximum value of h on P is attained at a unique point z = (Ci, • • • , Cn)- 
If X = (xi, . . . , x n ) is a vector of independent Bernoulli random variables such 
that Exj = Qj then the value of the probability mass function X at a point x G 
P fl {0, l} n is equal to exp {— hg(z) + £(x)}. In particular, 

^ exp{£(x)} = exp{h e (z)}Pr{X E P). 

xePn{o,i}" 

(3.4) Comparison with the Monte Carlo method. Suppose we want to sam- 
ple a random 0-1 point from the uniform distribution on Pfl {0, l} n . The standard 
Monte Carlo rejection method consists in sampling a random 0-1 point x, accept- 
ing x if x G P and sampling a new point if x ^ P. The probability of hitting P 
is, therefore, 2 _n |PflZ n |. It is easy to see that the largest possible value of h 
in Theorem 3.3 is nln2 and is attained at £i = . . . = £ n = 1/2. Therefore, the 
rejection sampling using the maximum entropy Bernoulli distribution of Theorem 
3.3 is at least as efficient as the standard Monte Carlo approach and is essentially 
more efficient if the value of h(z) is small. 

Applying a similar logic as in Section 3.2, we obtain the Gaussian heuristic 
approximation of (2.5.1). 

We notice that 

is the entropy of the Bernoulli distribution with expectation £ while 

<7(£) = (£+l)ln(£ + l)-£ln£ 

is the entropy of the geometric distribution with expectation £. One can suggest 
the following general maximum entropy approach, cf. also a similar computation 
in [Ja57]. 
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(3.5) Theorem. Let S C R n be a finite set and let conv(S') be the convex hull 
of S. Let us assume that conv(S') has a non-empty interior. For x G conv(.S'), let 
us define 4>{x) to be the maximum entropy of a probability distribution on S with 
expectation x, that is, 



<f>(x) = max 'S^Ps In — 

its Pi 

Subject to: ^^Ps = 1 



ses 

^2sp s = X 

ses 

Ps > for all s G S. 

Then <p(x) is a strictly concave continuous function on conv(S) . 

Let A C M, n be an affine subspace intersecting the interior of conv(S') . Then cj) 
attains its maximum value on A n conv(S') at a unique point z in the interior of 
conv(S'). There is a unique probability distribution p on S with entropy <f>(z) and 
expectation in A. Furthermore, the probability mass function of p is constant on 
the points of S fl A and equal to e~^ z ^ : 



In particular, 



Proof. Let 



^{s} = e ~^ z) for all s G S n A. 
\SnA\ = e^ z) p{Sn A}. 

h[p s : seS^J^Ps In — 

ses Ps 



be the entropy of the probability distribution {p s } on S. 

Continuity and strict concavity of <f> follows from continuity and strict concavity 
of H. Similarly, uniqueness of p follows from the strict concavity of H. 

Since 

9 -tf = ln!-l, 



dp s p. 



s 



which is finite for p s > and is equal to +oo for p s = (we consider the right 
derivative), we conclude that for the optimal distribution p we have p s > for all 
s. 

Suppose that A is defined by linear equations 



(ai,x) = Pi for i = 1, . . . ,d, 
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where G M n are vectors, ft G I are numbers and (•, ■) is the standard scalar 
product in M n . Thus the measure n is the solution to the following optimization 
problem: 



-> max 

Ps 



SjPs In — 

Subject to: = 1 

ses 

^2(ai, s)p s = pi for i = l,...,d 

ses 

p s >0 for all s e S. 

Writing the optimality conditions, we conclude that for some Ao, Ai, . . . , A^ we have 

d 

lnp s = A + ^2\i{a,i,s). 
i=i 

Therefore, 

Ps = exp |a + ^2\i{ai, s)| . 
In particular, for s G A we have 



i=i 



p s = exp | A + | ■ 

On the other hand, 

<P(z) =H(p s : 8 G tf) 
/ d 

= - I A + M a *» s ) 

ses V i=i 
d 

= — Ao — ^ 

i=l 

which completes the proof. □ 

Finally, we discuss a continuous version of the maximum entropy approach. 
We recall that x is an exponential random variable with expectation ( > if the 
density function ip of x is defined by 



(1/C)e" r / C for r>0 

for r < 0. 

17 



We have 



E x = ( and var x = ( 2 . 



The characteristic function of x is defined by 

E e lTX = for t el. 

1 - i(t 

(3.6) Theorem. Let P C W 1 be the intersection of an affine subspace in MJ 1 and a 
non-negative orthant W+. Suppose that P is bounded and has a non-empty interior. 
Then the strictly concave function 

n 

f(x) = n + ^ln^- for x = (£1, . . . , f n ) 

3 = 1 

attains its unique maximum on P at a point z = . . . , ( n ), where Q > for 
j = 1,... ,n. 

Suppose now that Xj are independent exponential random variables with expecta- 
tions (j for j = 1, ... ,n. Let X = (x\, . . . ,x n ). Then the density of X is constant 
on P and for every x G P is equal to e~^ z ^ . 

Proof. As in the proof of Theorem 3.1, we establish that Q > for j = 1, . . . , n. 
Consequently, the gradient of / at z must be orthogonal to the affine span of P. 
Assume that P is defined by a system of linear equations 

n 

^ ctijij = Pi for i = 1, . . . , d. 

3=1 

Then 

1 d 

— = ^2\ i a ij for j = l,...,n. 

^ i=i 

Therefore, for any x G P, x = (£i, . . . , £ n ), we have 

d 

i=i 



din 



In particular, substituting £j = (j, we obtain 



3 = 1 ^ 
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Therefore, the density of X at x G P is equal to 





exp \ - y -j- ) = e 

□ 

Again, X has the maximum entropy distribution among all distributions on W] 
subject to the constraint EIgF. 

A similar formula can be obtained for the exponential integral 

J e^) dx, 

where £ : MJ 1 — > R is a linear function, 

£(x) = 7i£i + . . . + 7 n ^ n for x = (£1, . . . , f„) . 
The integral may converge even if P is unbounded. We introduce 

n 

fi(x) = ra + J^lnfj + 7jCj. 
j'=i 

If £ is bounded from above on P and attains its maximum on P on a bounded 
face then the maximum of fg on P is attained at a unique point z = (Ci, • • • , Cn)- 
If X = (xi, . . . , x n ) is a vector of independent exponential random variables such 
that E xj = Q then the density of X at a point x E P is equal to 

exp{ +£(*)}. 

(3.7) The Gaussian heuristic for volumes. Below we provide an informal 
justification of the Gaussian approximation formula (2.1.1) 

Let P be a polytope and let x±, . . . , x n be the random variables as in Theorem 
3.6. Suppose that P is defined by a system Ax = 6, x > 0, where A = (ctij) is a 
d x n matrix of rank d < n. Let Y = AX, so Y = (yi, . . . , yd), where 

n 

Hi = ^2 ai i x i for » = 1> • • • > 

In view of Theorem 3.6, the density of Y at b is equal to 

(volP)e-^) (det AA T )~ 1/2 

(we measure volP as the (n — d) -dimensional volume with respect to the Euclidean 
structure induced from M. n ). 

We have Ey = b. The covariance matrix Q = (g^) of 1" is computed as follows: 

n n 

q i:j = cov (yi,yj) = a ik a jk var x k = J^a ifc ay fc Cfc- 
fc=i fc=i 
Assuming that the distribution of Y - at y = b is well approximated by the Gaussian 
distribution, we obtain formula (2.1.1) 
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4. Volumes of multi-index transportation polytopes 

We apply Theorem 2.2 to compute volumes of multi-index transportation poly- 
topes. We begin our discussion with ordinary (two-index) transportation polytopes. 
Although Theorem 2.2 does not imply the validity of the Gaussian approximation 
here, two-index polytopes provide a simple model case of computations that we 
later use in the case of a larger number of indices. 

(4.1) Transportation polytopes. 

For integers to, n > 1 let us choose positive numbers R = (ri,... ,r m ) and 
C = (ci, . . . , c n ) such that 

r t + . . . + r m = ci + . . . + c n = N 

and let us consider the polytope P = P{R, C) of all to x n non-negative matrices 
x = with the row sums ri, . . . ,r m and the column sums c\, . . . ,c n . As is 
known, P is a non-empty (to — l)(n — l)-dimensional polytope, also known as a 
transportation polytope, see, for example, [Y+84]. If to = n and R = C = (1, . . . , 1) 
then P is the polytope of n x n doubly stochastic matrices, also known as the 
Birkhoff polytope. We note that the row and column sums are not independent, 
since the total sum of all row sums is equal to the total sum of the column sums. 
We define the affine span of P by the following non-redundant system of linear 
equations: 

n 

=n for i = 1, . . . , m - 1 

m 

(4.1.1) ^2 = c i for 3 = • • • ' n ~ 1 and 

i=i 

E ^ = n. 

l<i<m 
l<j<n 

In other words, we prescribe the sums of the first to — 1 rows, the first n — 1 columns, 
and the total sum of the matrix entries. We observe that every column a of the 
matrix A of the system (4.1.1) contains at most 3 non-zero entries (necessarily equal 
to 1), so ||a|| < a/3. 

Let z = (dj) be a matrix, z G P, maximizing 

f(x) = mn+ E ln ^j- 

l<i<m 
l<j<n 



Theorem 2.2 associates with system (4.1.1) the following quadratic form q defined 

20 



on R m+n_1 : 

q(a; b;u)=^ ^ (J (a* + Pj + <^Y 



2 

l<i<m-l 
l<j<n-l 

m— 1 

2 



i=l 
n— 1 



2 

7 = 1 

+ 2^ mn ' 



Here 



a = («i, . . . ,a m _i) , 6 =(/?i,... ,/? n _i) 



are real vectors and a; is a real number, so (a; 6; u) is interpreted as a vector from 

j^m+n— 1 

To bound the eigenvalues of q from below, we bound the eigenvalues of a simpler 
form 



l<i<m-l 
l<j<n-l 

m—1 n — 1 

i=l j=l 

Let us consider the (m — 2)-dimensional subspace H a C ]R m + n_1 defined by the 
equations 

m— 1 



ctj = 0, 6 = 0, and to = 



and the (n — 2)-dimensional subspace H b C IR m+n 1 defined by the equations 



n-l 



/?j =0, a = 0, and u; = 0. 



i=i 



We observe that H a is an eigenspace of q with the eigenvalue n (since the gradient 
of q at x G H a is equal to 2nx) and that Hb is an eigenspace of q with the eigenvalue 
m (since the gradient of q at x G Hb is equal to 2mi). Let L C K m + n_1 be the 
orthogonal complement to H a + Hb. Then dimL = 3 and L consists of the vectors 

q» • • • ; « ; £» • • • ; ff ; ^ 

Vm-1 times n—1 times 
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for some real ct,f3, and u. Therefore, the restriction of q onto L can be written as 



a, .. .,«;/?,.. . ,/3;u 
„ ' > „ ' 

\m— 1 times n — 1 times 



(m - l)(ra - l)(a + + u) 2 + (m - l)(a + w) 2 + (ra - + w) 2 + w 2 . 

Since 

(a + /3 + w) 2 + (a + w) 2 + (/3 + w) 2 > 5 (a 2 + (3 2 + u 2 ) 

for some absolute constant 5 > and all a, /3 and a;, we conclude that the eigen- 
values of q exceed 

min{m — 1, n — 1} 
o - 

maxjm — 1, n — 1} 

for some absolute constant S > 0. Same holds for the eigenvalues of q as long as 
the numbers Qj are uniformly bounded away from 0. 

We notice that the minimum eigenvalue of q is too small to satisfy the conditions 
of Theorem 2.2. In fact, as Canfield and McKay have shown [CM07b], the volume 
of the Birkhoff polytope is not asymptotically Gaussian as m = n — > +oo, since 
there is a fourth-order correction akin to the Edgeworth correction. However, a 
very similar analysis can be applied to certain higher-dimensional versions of trans- 
portation polytopes and there it produces more satisfying results: asymptotically, 
volumes of such polytopes turn out to be given by the Gaussian formula (2.1.1). 

(4.2) Multi-index transportation polytopes. Let us fix an integer v > 2 and 
let us choose integers ki, . . . ,k u > 1. We consider the polytope of P of ki x . . . x k v 
arrays of no n- negative numbers where 1 < j% ^ ki for i = 1,... ,v, with 

prescribed sums along the coordinate hyperplanes. Namely, we choose positive 
numbers fy, where 1 < j < k^ for i = 1, . . . , v and such that 

for some iV and alH = 1, . . . , v and define P by the inequalities 

> for all ji,...,j v 

and equations 

jl,— ,3i-l,3i+l,— ,3v 

(4.2.1) for i = 1, . . . , v and I < j < h — 1 and 

jl, — ,jv 
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Let us choose a pair of indices 1 < i < v and 1 < j < hi — 1. We call the first sum in 
(4.2.1) the j-th sectional sum in direction i. Hence for each direction i = 1, . . . , v 
we prescribe all but the last one sectional sum and also prescribe the total sum 
of the entries of the array. When v = 2 we obtain the transportation polytope 
discussed in Section 4.1 We observe that every column a of the matrix A of the 
system (4.2.1) contains at most v + 1 non-zero entries (necessarily equal to 1), so 

INI < + 1. 

Let z = {(,j 1 ...j v ) be the point maximizing 

f{z) = kf-k v + 1,1 £./... x 

jl,--- ,jv 

on P. We describe the quadratic form q : M. d — ► R which Theorem 2.2 associates 
with system (4.2.1). We have d = k± + . . . + k v — v-\-l and it is convenient to think of 
R d as of a particular coordinate subspace of a bigger space V = R fcl © . . . ©R fc " ©R. 
Namely, we think of V as of the set of vectors (£, a;), where 

* = ( r ij) f° r 1 — J ' — ki an< ^ i = 1, • • • , v 

and and u> are real numbers. We identify R d with the coordinate subspace 
defined by the equations 

Tlfci = T2fc 2 = • • • = t v \- v =0. 
Next, we define a quadratic form p : V — > R by 

ii , • • • 

Then the quadratic form q of Theorem 2.2 is the restriction of p onto R d . 

To bound the eigenvalues of q from below, we consider a simpler quadratic form 
q which is the restriction of 

p(t, <jj) = ( r iJi + • • • + T U j v + uf 

jl,— ,jv 

onto R d . 

For % = 1, . . . , v, let us consider the (ki — 2) -dimensional subspace Hi C R d 
defined by the equations 

fci-i 

Tij = 0, Ti'j = for i' ^ i and all j, and a; = 0. 
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Then Hi is an eigenspace of q with the eigenvalue 



/ 



\ — fci • • • ki—iki-^-i ■ ■ ■ kjy, 

since the gradient of q at x e Hi is equal to 2XiX. Let L C M d be the orthogonal 
complement to Hi © . . . © H v in M d . Then dimL = u + 1 and L consists of the 
vectors 

/ \ 

on, ■ ■ ■ , Q-i , 0; «2, ■ • • , «2 , 0; . . . , o^, . . . , ay ; 0; u 

Y fci — 1 times &2 — 1 times fc„ — 1 times 

for some real a±, . . . , a v ; u>. Denoting 

A*o = (ki — 1) • • • (A;^ — 1) and 

^ = (fci - 1) • • • (ki-! - 1) (fc i+ i - 1) • • • (k v - 1) 

We observe that the restriction of q onto L satisfies 
/ \ 



Q 



ail, • • • ,cti , 0; a.2, ■ • • , Q-2 , 0; . . . , ay, . . . , ay ; 0; cj 

y fci — 1 times &2 — 1 times k v — l times y 



2 ^ — ^ r 

> («! + • • • + OL v + <Jj) + A*i ("1 + • • • + + "i+l + • • • + a u + U Y 



1=1 



Since 



i=l 



2 ^ — r 

(ai + . . . + a u + u) + (ai + ... + a.i-i + a i+ i + . . . + a u + u) 
> 6 ( + 



i=i 



for some S = 8(v) > and all oti, . . . , ay and u, we conclude that the eigenvalues 
of q exceed 



S(y) min (ki - 1) 2 FT (kj - 1) , 
i=i,...,f J-- 1 - 



where 5(z/) > is a constant depending on z/ alone. 

Suppose now that v is fixed and let us consider a sequence of polytopes P n 
where fci, . . . , Ay grow roughly proportionately with n and where the coordinates 
Cji—jv remain in the interval between two positive constants. Then the minimum 
eigenvalue of the quadratic form q in Theorem 2.2 grows as O (n^~ 2 ) . In particular, 
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for v > 5 Theorem 2.2 implies that the Gaussian formula (2.1.1) approximates the 
volume of P n with a relative error which approaches as n grows. 

As an example, let us consider the (dilated) polytope Pk of polystochastic tensors, 
that is k x ... x k arrays of non-negative numbers with all sums along coordinate 
hyperplanes equal to k v ~ x , cf. [Gr92]. By symmetry, we must have 

Cji...j„ = i- 

Theorem 2.2 implies that for v > 5 

vo\P k = (l + o(l)) (2?r)(; ^ +1)/2 as £;^+oo. 

Interestingly, for v = 2, where our analysis is not applicable, the formula is smaller 
by a factor of e 1//3 than the true asymptotic value computed in [CM07b]. 

5. The number of multi-way contingency tables 

We apply Theorems 2.4 and 2.6 to compute the number of multi-way contingency 
tables. The smallest eigenvalue of the quadratic form q is bounded as in Section 
4 and hence our main goal is to bound the additive error A. Again, we begin our 
discussion with ordinary (two-way) contingency tables, where Theorems 2.4 and 2.6 
do not guarantee the validity of the Gaussian approximation, but which provide a 
simple model case for computations used later in the case of multi-way tables. 

(5.1) Contingency tables. Let us consider the transportation polytope P = 
P{R,C), see Section 4.1, where the row sums r 1? ... , r m and the column sums 
ci, . . . , c n are integer. Integer points in P(R, C) are called contingency tables and 
0-1 points in P(R, C) are called binary contingency tables with margins R and C, 
see [DE85]. 

We assume that P is defined by system (4.1.1). To estimate the additive error 
term A in Theorems 2.4 and 2.6, we need to construct sets of integer vectors of the 
following three types: 

for k = 1, . . . , m — 1 we construct a set of m x n integer matrices y such that 
the k-th row sum of y is 1, all other row and column sums, with possible exceptions 
of the m-th row sum and n-th column sums are 0, and the total sum of the matrix 
entries is as well; 

for k = 1, . . . , n — 1 we construct a set Y^ of m x n integer matrices y such 
that the k-th column of y sum is 1, all other row and column sums, with possible 
exceptions of the m-th row sum and the n-th column sums are 0, and the total sum 
of the matrix entries is as well; 

we construct a set Y$ of m x n integer matrices y such that all the row and column 
sums of y with possible exceptions of the m-th row sum and the n-th column sum 
are 0, and the total sum of the matrix entries is 1. 
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To construct Y^, let us choose an integer 1 < I < n and let us define a matrix 
V = (Vij) by letting r/ki = 1, f]mi = — 1 and letting all other entries rjij equal to 
0. The set Y^ 1 contains n vectors y with pairwise disjoint support and hence the 
maximum eigenvalue p^} of the corresponding quadratic form 



1 n 

~ Yl (Zki -£,mi) 2 for X = (Cij) 



n 
i=i 



is 2/n. Similarly, to construct Y£ , let us choose an integer 1 < I < m and let us 
define a matrix y = (rjij) by letting rjik = 1, r\\ n = — 1 and letting all other entries 
rjij equal to 0. The maximum eigenvalue p% of the corresponding quadratic form 



vf(*) = -X>*-fc») s 



m 
i=i 

is 2/m. 

Finally, to construct Yq, let us choose two indices 1 < k < m — 1 and 1 < I < n — 1 
and let us define a matrix y = (rjij) by letting rjki = —1, rjk n = 1 5 Vmi = 1 and 
letting all other entries rjij equal to 0. For the corresponding quadratic form tjj we 
have 



(m - l)(n - 1) f-r 1 

l<fc<m — 1 
l</<n-l 

-( m -l)(n-l) E 3 (&» + ^mi + ti) 

v ;v ; Kfe<m-1 



KKn-1 
Q m — 1 _ n— 1 



m-lS^ + n _iE^+( m _i)( n _i) E 

fc=l /=1 v ;v ' Kfc<m-1 



KKn-1 



Hence we can choose 

/9 = max 



m — 1 ' n — 1 



in Theorems 2.4 and 2.6, so the additive term A is exponentially small in min{m, n}. 
This bound is pretty weak but it is getting better as we pass to multiway tables. 
In fact, as Canfield and McKay have shown [CM07a], in the simplest case of R = 
(r, . . . , r) and C = (c, . . . , c), the number of contingency tables is not given by the 
Gaussian formula, since there is a 4-th order term correction. 
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(5.2) Multi-way contingency tables. Let us consider the z/-index transporta- 
tion polytope P of Section 4.2. We assume that the affine span of P is defined 
by system (4.2.1), where numbers fiij are all integer. The integer points in P are 
called sometimes multi-way contingency tables while 0-1 points are called binary 
multi-way contingency tables, see [Go63] and [DL05]. 

To bound the additive error term A in Theorems 2.4 and 2.6, we construct a 
set Yij of k\ x . . . x k v arrays y of integers such that the total sum of entries of y 
is 0, the j-th sectional sum in the i-th direction is 1 all other sectional sums are 
0, where by "all other" we mean all but the /cj-th sectional sums in every direction 
% = 1, . . . , v. For that, let us choose v — 1 integers mi, . . . , mj_i, mj+i, . . . , m v , 
where 

1 < mi < ki, . . . , 1 < rrii-i < ki-i, 1 < m^i < fc^+i, . . . , 1 < m v < k v 
and define y = (r]j 1 ...j l/ ) by letting 

Vm 1 ...m i - 1 ,j,m i+1 ...m v = 1) Vm 1 ...mi- 1 ,ki,m i+1 ...m v = ~ 1 

and letting all other coordinates of y equal to 0. 

Thus the set contains k\ • • • ki-iki+i • • • k v elements y, and the corresponding 
quadratic form tpij can be written as 

V'ij ( x ) = |v I {£m 1 ---m i - 1 ,j,m i+1 ---m„ ~ £m 1 ---m i - 1 ,k i ,m i+1 ---m„) 

\ i 7 

■> mi,...,mj-i,mj|i,...,m„ 

for x= 

from which the maximum eigenvalue /9y of is 2/ki ■ ■ ■ ki-iki + i • • ■ k v . 

Next, we construct a set Y of arrays y of k\ • • • k v integers such that 

the total sum of entries of y is 1 while all sectional sums with a possible exception 
of the ki-th sectional sum in every direction i are equal 0. For that, let us choose 
v integers mi, . . . , m^, where 

1 < mi < k\ — 1, . . . , 1 < m v < k v — 1 

and define y = (// /: ) by letting 

Vm\...ra v 1 ^ 
yki,m 2 -..m l , = 1 
ym 1 ,k2,m 3 ...m v = 1 

ymi...m v -i,k u = 1 



and by letting all other coordinates equal to 0. 
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The set Y contains (ki — 1) ■ • ■ (k v — 1) elements and the corresponding qua- 
dratic form ip of Theorems 2.4 and 2.6 can be written as 



y • ^ 



1 \ 2 

/ j ((1 ~~ l/ )£mi...m„ + 6i,m2...ni„ + . . . + £ mi . . .m„„ ! ,fc„ ) 



\ Y 0, 

l<7Tll<fcl— 1 

\<m v <k v — 1 



— iy I / , , ^v- 1 - w Smi-.m,, T Si; 1 ,m2...m„ T • • ' T Smi...m„_i,l;„^ 

Kmi<fei- 1 

l<m„ <fc„ — 1 

Therefore, the maximum eigenvalue po of V'o does not exceed 

(z/ + — l) 2 max <— r — — r 

and the same bound can be used for the value of p in Theorems 2.4 and 2.6. 

Suppose now that v is fixed and let us consider a sequence of polytopes P n where 
fei, . . . , &v grow roughly proportionately with n. Then in Theorems 2.4 and 2.6 we 
have 

= o< 1 



n v-\ 



Let us apply Theorem 2.6 for counting multi-way binary contingency tables. We 
assume, additionally, that for the point z = (Cji...j„) maximizing 

/(*>= £ j J - ln r- + r^— 

on the transportation polytope _P n we have 

— >ji i • • • — 

for some constant 1/2 > 5 > and all ji, . . . , .?'„• Then we can bound the additive 
term by 

|A| < expl-7^- 1 } 
for some constant 7 > 0. On the other hand, by Hadamard's inequality, 

det BB T = n 0{n K 

Therefore, for v > 3, the additive term A is negligible compared to the Gaussian 
term. From Section 4.2, we conclude that for v > 5 the relative error for the num- 
ber of multi-way binary contingency tables in P n for the Gaussian approximation 
formula (2.5.1) approaches as n grows. 
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Similarly, we apply Theorem 2.4 for counting multi-way contingency tables. Here 
we assume, additionally, that for the point z = maximizing 

f( X ) = Yl foi-J* + X ) ln fai-Jv £/:•.../• ..,/■ 

3i,- ,Jk 

on the transportation poly tope _P n the numbers Q 1 ...j v lie between two positive 
constants. As in the case of binary tables, we conclude that for v > 3, the additive 
error term A is negligible compared to the Gaussian approximation term as n — > 
+oo. Therefore, for v > 5 the relative error for the number of multi-way contingency 
tables in P n for the Gaussian approximation formula (2.3.1) approaches as n 
grows. 

Computations show that in the case of k\ = . . . = k v = k for the matrix A of 
constraints in Theorems 2.4 and 2.6 we have 

det AA T = jfc(" 2 -")(fc-D. 

Hence we obtain, for example, that the number of non-negative integer z/-way k x 
. . . x k contingency tables with all sectional sums equal to r = ak v ~ x is 

(l + o(l)) {(a + l) a+1 a-«f (27ra 2 + 27ra)- (fc ^ +1)/2 ^— 2 )( fc - 1 )/ 2 

provided v > 5, k — > +oo and a stays between two positive constants. Interest- 
ingly, for v = 2 (where our analysis is not applicable) the obtained number is off 
by a constant factor from the true asymptotic obtained in [CM07a] . 

Similarly, the number of binary z/-way k x . . . x k binary contingency tables with 
all sectional sums equal to r = ak v ~ l is 

(l + o(l)) (^(l-a) 1 -")"^ (2rra-27ra 2 )- (fc ^ +1)/2 ^ 2 )( fc - 1 )/ 2 

as long as v > 5, k — ► +oo and a remains separated from and 1. Again, for v = 2 
the formula is off by a constant factor from the asymptotic obtained in [C+08]. 

6. Proof of Theorem 2.2 

We start with some standard technical results. 

(6.1) Lemma. Let x±, . . . ,x n be independent exponential random variables such 
that ~EiXj = Cj for j = 1, . . . ,n, let ai, . . . , a n E M. d be vectors which span M. d and 
let Y = x\a\ + . . . + x n a n . Then the density ofY at b £ is equal to 

_J_ f e -*{b,t) I fr i 

(27ryJ Rd [j^l-zQia^t} 
Proof. The characteristic function of Y is 

Ee ^) = ff J- -. 

The proof now follows by the inverse Fourier transform formula. □ 
We need some standard estimates. 
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(6.2) Lemma. Let q : R d — > R be a positive definite quadratic form and let u > 
be a number. 

(1) Suppose that ui > 3. Then 

[ e-«(*) eft < e"^ 2 / e-«W dt. 

it: q(t)>uid JR d 

(2) Suppose that for some X> we have 

q(t) > X\\t\\ 2 for all t e R d . 
Let a e R d be a vector. Then 

e- q W dt < e~ x ^ I e~ q ^ dt. 



't: \{a,t)\>u\\a\\ 

Proof. We use the Laplace transform method. For every 1 > a > we have 



/ e q( - ' dt < exp {a(q(t) -ud) -q(t)} dt 

Jt: a(t)>ujd Jt: a(t)>Lod 

/ exp{-(l -a)q(t)} dt 
e- q V dt 



•t: q(t)>uid Jt: q(t)>u>d 

< e~ aujd 

„—aujd 



(l-a) d / 2 

Optimizing on a, we choose a = 1 — l/2u to conclude that 

/ e~ q ^ dt < exp(-wd+^ + ^hn>)l / e'^ dt. 

Jt: q(t)>u;d { 2 2 J J R d 

Since 

ln(2w) < u - 1 for u > 3, 

Part (1) follows. 

Without loss of generality we assume that a ^ in Part (2). Let us consider the 
Gaussian probability distribution on R d with the density proportional to e~ q . Then 
z = (a,t) is a Gaussian random variable such that Ez = and varz < ||a|| 2 /2A. 
Part (2) now follows from the inequality 

Pr {\y\ > t) < e" r2/2 

for the standard Gaussian random variable y. □ 
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(6.3) Lemma. For p > and k > d we have 

/ (1 + 11*11) dt < (1 + p ) 

JteR d : \\t\\> P 1 - d ) 

Proof. Let C M d be the unit sphere in M, d . We recall the formula for the 

surface of S"" 1 : 

1 1 r(d/2) 

We have 

r> + CSO 

f||>/' "V 



(l + ||t|| 2 )- fc/2 dt = \§ d ~ 1 \ / (l + r 2 )-" 72 ^" 1 dr 

„£K d : \\t\\>p Jp 

<\^\ / + °°(l + r 2 ) (d - fc - 2)/2 rdr, 
Jp 

where we used that 

r d-l =rr d-2< r(l + r 2 ) (rf-2)/2_ 

The proof now follows. □ 

Now we are ready to prove Theorem 2.2. 

(6.4) Proof of Theorem 2.2. Scaling vectors cij if necessary, without loss of 
generality we may assume that 9 = 1. 

From Section 3.7 and Lemma 6.1, we have 

volP = e^(detAA T ) 1/2 -i- / e -*<*,*> [ TT J_ ) dt 

Hence our goal is to estimate the integral and, in particular, to compare it with 

/ dt = (27v) d / 2 {detBB T )~ 1/2 . 

Let us denote 



for t e 



Let 

a = 4d+ 10 In-. 

e 
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We estimate the integral separately over the three regions: 



the outer region ||£|| > 1/2 

the inner region q(t) < a 

the middle region ||t|| < 1/2 and q(t) > a. 

We note that for a sufficiently large constant 7 we have q(t) > a in the outer 
region, we have \\t\\ < 1/2 in the inner region and the three regions form a partition 
ofR d . 

We start with the outer region > 1/2. Our goal is to show that the integral is 
negligible there. 
We have 

1/2 



Let us denote 

tj = (j( a j, t ) 2 for j = l,.-.,n. 
The minimum value of the log-concave function 

n 

on the polytope 

n 

> 2 M\t\\ 2 an d < €3 < \\tf 

i=i 

is attained at an extreme point of the polytope, that is, at a point where all but 
possibly one coordinate £j is either or ||£|| 2 . Therefore, 



1/2 

2 A+l/2 



t) 2 

Applying Lemma 6.3, we conclude that 



< (l + ll^ll 2 )" 



ha*-. imi>i/2 1 {t)l 1 ~ r(d/2)(2A-d-i) U; 



27r d/2 /g\ (d-2A+l)/2 

l>i/s 

By the Binet-Cauchy formula and the Hadamard bound, 



det BB T < Q ] < n" 
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It follows then that for a sufficiently large absolute constant 7 and the value of the 
integral over the outer region does not exceed (e/f 0)(27r) d / 2 det(SS T )- 1 / 2 . 

Next, we estimate the integral over the middle region with < 1/2 and q(t) > 
a. Again, our goal is to show that the integral is negligible. 

From the estimate 



in(i+o-e+ 



e e 



i£i 4 1 

< — - for all complex |£| < -, 



we can write 



ln(l -iQiaj.t)) = -iQ( aj ,t) + l^(aj,t) 2 + ^(a^t)* + gj {t)$( aj ,t)\ 



where 



Since 



we have 



(6.4.1) 



In particular, 



\9j(t)\ < - for j = 1,... ,n. 



F(t)=exp{-q(t)-if(t) + g(t)} 

1 n 

where /(f) = - £ $ (a, , t) 3 and 

i=i 



1 n 



\F(t)\ < e - 3q W/ 4 provided < 1/2. 
Therefore, by Part (1) of Lemma 6.2 we have 



F(t) dt 



ll*ll<l/2 
q(t)>a 



< j e -3«(t)/4 dt 

't: q(t)>a 



<e -3d/2 e 3 / e -3q(t)/4 ^ 



<e° / e~ q{t) dt. 
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Finally, we estimate the integral over the inner region where q(t) < a and, 
necessarily, < 1/2. Here our goal is to show that the integral is very close to 

e- q ^ dt. 
From (6.4.1), we obtain 



(6.4.2) 



F(t) dt- e~ q W dt 

t: q(t)<a Jt: q(t)<a 



< f e -9(t) 


e -if(t)+g(t) _ 1 


Jt: q(t)<a 




a then \\t 2 < a/A and hence 




i n 


n 



dt. 



3 = 1 



2A 



j=i 



A 



Thus for all sufficiently large 7, we have \g(t) \ < e/10. 
Let 

X = jt : q(t)<a and 01(^)1 < for j = l,...,n| 

By Part (2) of Lemma 6.2, for all sufficiently large 7, we have 



'R d \X 

whereas for t £ X we have 



f e -m dt < ± / e-" [t) dt 

JRd\ y 10 



n n 

1/(01 < «E<?i< a i»*)i 3 ^ ^E^*) 2 



3^^ IXJ,/I "30a^^ NJ '' -15 
3=1 3=1 



Estimating 

e -if(t)+g(t) _ 1 



< - for t e X and 



e -if(t)+g(t) _ x 



< 3 for t £ X 



we deduce from (6.4.2) that 



/ F(t) dt - [ e~ q{t) dt 

Jt: q(t)<a Jt: q(t)<a 



Since by Part (1) of Lemma 6.2, we have 



< 



3 / e-*W dt+\ ( e" 9 « dt 
Jm. d \x 3 J x 



< 



2e 



e~ q{t) dt. 



L 



e~ q ^ dt < e" 2d e 5 / e"«W ^ 



q(t)><r 



the proof follows. 



□ 
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7. Proof of Theorem 2.6 
First, we represent the number of 0-1 points as an integral. 

(7.1) Lemma. Letpj, qj be positive numbers such thatpj+qj = 1 for j = 1, . . . ,n 
and let p be the Bernoulli measure on the set {0, l} n of 0-1 vectors: 

n 

V>{ x } = II P J 'V f° r X = (&> • • • ' €n) ■ 

3 = 1 

Let P C M, n be a polyhedron defined by a vector equation 

£iai + . . . + i n an = b 
for some integer vectors a±, . . . , a n ; b E Z d and inequalities 

< ,£n < i. 

Let IT C W 1 be the parallelepiped consisting of the points t = (n, . . . , To) such that 

< T~k < tv for k = 1, . . . , d. 

Then, for 

up) = ^ 

x6Pn{o,i}™ 

we have 

i c n 

mp) = j^r L e ~ l(t,b) n (/''-'/''''"'• , ) dt 

Here (•, •) is the standard scalar product in M, d and dt is the Lebesgue measure on 



Proof. The result follows from the expansion 

n 'i 

II (Pj + Qje*^) = E exp{z<£ia 1 + ... + £nan, I}}]]!') S 



n 

3 

3=1 ' x6{0,l}~ 3 = 1 

X=(€l,--- ,£n) 



and the identity 

1 fun f 1 if u = 

(27r) d J n " 10 if«GZ d \ {0}. 



□ 

The integrand 

n 
3 = 1 

is the characteristic function of Y = AX where X is the multivariate Bernoulli 
random variable and A is the matrix with the columns a±, . . . ,a n . 
The following result is crucial for bounding the additive error A. 
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(7.2) Lemma. Let A be a d x n integer matrix with the columns ai, . . . , a n G Z d . 
For k = 1, . . . , d let Y\~ C Z n be a non-empty finite set such that Ay = e^ for all 
y G Ifc; where e^ is the k-th standard basis vector. Letifjk : R n — ► R be a quadratic 
form, 

ipk(x) = -r—r (y> x ) 2 f° r x e Rn ' 

yeY k 

and let pk be the maximum eigenvalue of of ifjk . 

Suppose further that < £i, . . . , ( n < 1 are numbers such that 

Q(l-Q)>a for some < a < 1/4. 

Then for t = (n, . . . ,Td) where —it < < 7r for k = 1, . . . , d we have 

n 

no <* • '"'•') 

i=i 

Proof. Let us denote 

n 

^)=n( i -o+c J ^^). 

Then 

n 

W)| 2 = II (C 1 - 0) 2 + 2C,(1 - C) cos(a„ f) + C 2 3 ) ■ 

3 = 1 

For real numbers £,77, we write 

£ = 77 mod 2-7T 
if £ — 77 is an integer multiple of 27r. Let 

— 7r < 7j < 7T for j = 1, . . . ,71 

be numbers such that 

(7.2.1) ( a ji~t) = 7j mod 2/t for j = 1, . . . , n. 

Hence we can write 

n 

I 2 = II (C 1 - C;) 2 + 2 0(! - Ci) cos 7, + C 2 ) • 

i=i 
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£ exp {-fi}' 



Since 

we have 
(7.2.2) 

Let 



7 

COS 7 < 1 — for — 7T < 7 < 7T, 

5 



< ft (l - 3^1?) <*J- 2 i±, 

j=l V ' / [ j=l 



c= (71,... , 7n ), c G W 1 . 

Then for all y G YJc we have 

Tk = (e*., £) = (Ay, t) = (y, = (y, c) mod 2tt, 
where A* is the transpose matrix of A. Since jr^l < 7r, we have 

|(y,c)| > |r fc | for all y G Yfe. 

Therefore, 



Icll 2 > 



c) 2 > ^. 



Pfc PfcKfel Pk 

The proof follows by (7.2.2). 

(7.3) Proof of Theorem 2.6. By Theorem 3.3 and Lemma 7.1, we write 

e h(z) 



□ 



(7.3.i) \p n {o, in = y n e ~ t{b,t) Hi 1 - o + v,<"' •' ) 



where IT is the parallelepiped consisting of the points t = (ti, . . . , r^) with — tv < 
Tk < 7r for A; = 1, . . . , d. 
Let us denote 



F(t) = e" l(M) n (l - Cj + Cje l(aj ' t} ) ■ 



3 = 1 



If 



we have 



< — , 
- 40' 



(a j5 t)|<- for j = l,...,n. 



Using the estimate 



< — for all real £, 
- 24 
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we can write 



e ^^) = l + i(a,,t)-i^-zi^+^(t)(a,,t) 4 , 



where \g 3 {t) I < ^ for j = 1,... ,n. 



Therefore, 



F(t) = e -«M> n (l + <>^> " O^f^ - iCi^^ + Qg^a^A . 

7 = 1 ^ ' 



J 

Furthermore, using the estimates 



£2 tS 

ln(l+0-e+y-y 



|4 



< i^J_ for all complex |f | < 1/2 



and that 



we can write 

F(t) = e~ q{t)+if{t)+9{t \ 
1 n 

where /(f) = - ^(2^ - 1) (Q - C?) (a,-, t) 3 and 
(7.3.2) J=1 

n 

b(t)|<2^(a„t) 4 . 
i=i 

In particular, 

|</(*)| < ^(t) provided < -L 

Let 

a = 4d+ 10 In-. 

e 

We split the integral (7.3.1) over three regions. 
The outer region: 

ll*l|oo>^. 

We let 



A =(2^F / F(i) *' 

ten 

l|t||oo>l/40 
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and use Lemma 7.2 to bound |A|. 
The middle region: 

q(t) > a and < 

From (7.3.2) we obtain 

\F(t)\ <e" 3 ^)/ 4 

and as in the proof of Theorem 2.2 (see Section 6.4), we show that the integral 
the region is asymptotically negligible for all sufficiently large 7. 
The inner region: 

q(t) < a. 

Here we have 

M. < Ml < 

provided 7 is sufficiently large. 



If q(t) < a then ||t||oo < 11*11 < a/ct/A and 

|^)| < 2$>,t> 4 < 2^^(0-C|)(a„t) 2 <4^. 
i=i j=i 

In particular, if constant 7 is large enough, we have \g(t) \ < e/IO. 
As in Section 6.4, we define 

X={t: q{t)<a and \{a jt t)\ < ^- for j = f,...,n| 

Hence for t e X we have 

1 n n 

< 6Eto-Ci)K«i»*>i 3 < 6^Eto-Ci)<«i»*> 2 < , - 

By Part (2) of Lemma 6.2, for all sufficiently large 7, we have 

e -9(t) d t < — [ e ~i(t) dt 

and the proof is finished as in Section 6.4. 

8. Proof of Theorem 2.4 

We begin with an integral representation for the number of integer points. 
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e 



(8.1) Lemma. Letpj, qj be positive numbers such thatpj+qj = 1 for j = 1, . . . ,n 

and let \x be the geometric measure on the set Z™ of non-negative integer vectors: 

n 

^i x ) = n Pjdj* x = (&» • • • ' &*) • 

3=1 

Let P C W 1 be a polyhedron defined by a vector equation 

£iai + . . . + i n a n = b 
for some integer vectors a±, . . . , a n ; b G Z d and inequalities 

Ci,... ,en>o. 

Le£ II C M d 6e i/ie parallelepiped consisting of the points t = (n, . . . , r^) swc/i t/iat 

— 7T < r/c < 7T /or fc=l,...,d. 

Then, for 

p(p)= J2 ti*} 

we have 

ffere (•, •) is the standard scalar product in M, d and dt is the Lebesgue measure in 
R d . 

Proof. As in the proof of Lemma 7.1, the result follows from the multiple geometric 
expansion 

n n 

IT i n P Ma„t) = E exp{i(eiai + ... + e„a„, *>}II^- 
z=(£i,--- 



□ 



The integrand 



Pj 



n i _ flf, e *< a i>*> 



is, of course, the characteristic function of Y = AX, where X is the multivariate 
geometric random variable and A is the matrix with the columns a\, . . . , a n . 
The following result is an analogue of Lemma 7.2. 
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(8.2) Lemma. Let A be a d x n integer matrix with the columns ai, . . . , a n E 7h d . 
For k = 1, . . . , d let Yj~ C Z d be a non-empty finite set such that Ay = for all 
y G Yfc, where is the k-th standard basis vector in Z d . Let tpk : M n — > M. be a 
quadratic form, 



and £e£ 6e £/ie maximum eigenvalue of tpk. Suppose further that Ci> • • • Xn > 
are numbers such that 

Q(l + Q) > a for some a > 0. 
Then for t = (n, . . . , r^) where —tv < Tk < tv for k = 1, . . . , d, we have 



1 



< 1 + -an' 



-m k 



where mu = 



Proof. Let us denote 



n 1 

= n i + o - 



Then 



Let 



1 

l F ( t )| 2 = n i + 2 Cj (l + Cj )(l-cos(a J ,t))' 
— 7r < 7j < 7T for j = 1, . . . , n 



be numbers such that 



7j = (cij, t) mod 27r for j = 1, . . . , n. 



Hence we can write 



1 

|F(t) ' 2 = 5 l + 2C J (l + C)(l-cos 7j ) 



< 



1 

11 1 + 2a(l - 



^ 1 + 2a(l — cos7j) ' 



Since 



7 

COS 7 < 1 for — 7T < 7 < 7T, 

5 
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we estimate 



(8.2.1) \F(t)\ 2 < n(l + jj«7j) 1 

Let 

c= (li,-- - ,7n)- 
As in the proof of Lemma 7.2, we obtain 



2 



INI 2 > -• 
pk 

Let us denote £j = 7? for j = 1, . . . , n. The minimum of the log-concave function 
on the polytope defined by the inequalities < £j < 7r 2 for j = 1, . . . , n and 



n _ 2 



is attained at an extreme point of the polytope, where all but possibly one coordi- 
nate £j is either or it 2 . The number of non-zero coordinates £j is at least r 2 /pk^ 2 
and the proof follows by (8.2.1). □ 

(8.3) Proof of Theorem 2.4. By Theorem 3.1 and Lemma 8.1, we have 



(8.3.1) |PnZ n | = ^ 7 [ e'^f] 7 - 



where IT is the parallelepiped consisting of the points t = (ti, . . . , r^) with — 7r < 
Tfc < 7r for /c = 1, . . . , d. 
Let us denote 



1 

Fm = e~ i{t ' b) TT - 



Similarly to the proof of Theorem 2.6 (see Section 7.3), assuming that ||t||oo < 1/46*, 
we write 



_ e -«(i)-i/(t)+9(t) ) 

1 n 

where /(f) = - £ + C*) (2£ + l)(a jt tf and 



\g(t)\ < 2^(1 + C,) 4 (a,^) 4 - 
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We let 

a = 4d + 10 In - 

e 

and as in the proof of Theorem 2.6 (see Section 7.3), we split the integral (8.3.1) 
over the three regions: 

the outer region: ||i||oo > 1/46*, 

the middle region: q(t) > a and ||£||oo < 1/40 and 

the inner region: q(t) < a. 



For the outer region, we let 



1 

(27) 



ten 
ll*lloo>i/4e 



and use Lemma 8.2 to bound A. 
We have 

\F(t)\ < e - 3 «(*)/ 4 

in the middle region and we bound the integral there as in Section 7.3. 
In the inner region, we have ||£||oo — \\t\\ < a/ct/A and 

\g{t)\ < 2^(1 + C i ) 4 (a„t) 4 < 2^^(C J+ C|)(a„t) 2 < 4°^-. 

3=1 3=1 

We define 

X = {t: q{t)<a and (2Q + 1) \( aj , t)\ < ^- for j = l,...,n} 
and note that for Jslwe have 

1 n n 
3=1 3=1 

The proof is finished as in Section 7.3. □ 
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